#include "ostu/ostu.h"

ostu::ostu(std::string imagePath){
	originImage = cv::imread(imagePath, 0);
	std::cout<<"size:"<<originImage.size<<std::endl;
	std::cout<<"type:"<<originImage.type()<<std::endl;
}


void ostu::drawHistogram(cv::Mat image){
	const cv::Mat originImageCopy = image.clone();
/*	auto countPixel = cv::Mat::zeros(1, 256, CV_32SC1);
	auto pixel = cv::Mat(1, 256, CV_8SC1);
	for(int j = 0; j < pixel.cols; j++){
		pixel.at<short>(0,j) = j;
	}*/
	std::vector<int> countPixel(256, 0);
	//直方图灰度分布概率
	std::vector<float> pixelProb;
	std::vector<int> pixel(256);
	for(int i = 0; i < 256; i++){
		pixel[i] = i;
	}


	for(int i = 0; i < originImageCopy.rows; i++){
		for(int j = 0; j < originImageCopy.cols; j++){
			countPixel[originImageCopy.at<uchar>(i, j)]++;
		}
	}

	float sum_pixel = originImageCopy.rows * originImageCopy.cols;
	for(auto iter = countPixel.begin(); iter != countPixel.end(); iter++){
		pixelProb.push_back(*iter /sum_pixel);
	}
	//plt::bar(pixel, countPixel);
	plt::stem(pixel, pixelProb);
	plt::show();
}


cv::Mat ostu::equalizeHistogram(){
	float sumPixel = originImage.rows * originImage.cols;

	int size = 256;

	auto tempImage = originImage.clone();
	//统计像素值出现的次数
	std::vector<int> countPixel(256, 0);
	//直方图灰度分布概率
	std::vector<float> pixelProb;
	//直方图均衡变换后得到的值
	std::vector<int> transformValue;
	//直方图均衡结果
	std::vector<int> equalizeResult(256, 0);

	//统计灰度值出现次数
	for(int i = 0; i < originImage.rows; i++){
		uchar* data = originImage.ptr<uchar>(i);
		for(int j = 0; j < originImage.cols; j++){
			countPixel[data[j]]++;
		}
	}

	//计算灰度值概率
	for(auto iter = countPixel.begin(); iter != countPixel.end(); iter++){
		pixelProb.push_back(*iter / sumPixel);
	}

	//计算直方图均衡结果
	for(int i = 0; i < pixelProb.size(); i++){
		int s = 0;
		for(int j = 0; j <= i; j++){
			s += round(255 * pixelProb[j]);
		}
		if(s > 255)
			s = 255;
		transformValue.push_back(s);
	}

	//进行直方图均衡变换
	for(int i = 0; i < tempImage.rows; i++){
		for(int j = 0; j < tempImage.cols; j++){
			tempImage.at<uchar>(i, j) = transformValue[tempImage.at<uchar>(i, j)];
		}
	}

	drawHistogram(tempImage);

	return tempImage;
}


void ostu::opencvOSTUtest(){
	auto tempImage = originImage.clone();
	cv::threshold(tempImage, tempImage, 0, 255, CV_THRESH_OTSU);
	cv::imshow("OSTU", tempImage);
	cv::waitKey(0);
}


void ostu::ostuAchieve(){
	auto tempImage = originImage.clone();
	auto tempImage_2 = originImage.clone();
	std::vector<int> countPixel(256, 0);
	//直方图灰度分布概率
	std::vector<float> pixelProb;
	std::vector<int> pixel(256);
	for(int i = 0; i < 256; i++){
		pixel[i] = i;
	}


	for(int i = 0; i < tempImage.rows; i++){
		for(int j = 0; j < tempImage.cols; j++){
			countPixel[tempImage.at<uchar>(i, j)]++;
		}
	}

	float sum_pixel = tempImage.rows * tempImage.cols;
	for(auto iter = countPixel.begin(); iter != countPixel.end(); iter++){
		pixelProb.push_back(*iter / sum_pixel);
	}

	//类内方差
	std::vector<float> class_within_variance_vec;
	//类间方差
	std::vector<float> class_between_variance_vec;

	for(int k = 0; k < 256; k++){
		//类0的概率
		float class_0_prob = 0;
		//类1的概率
		float class_1_prob = 0;
		for(int i = 0; i <= k; i++){
			class_0_prob += pixelProb[i];
		}
		class_1_prob = 1 - class_0_prob;

		//std::cout<<"class_0_prob:"<<class_0_prob<<" class_1_prob:"<<class_1_prob<<std::endl;

		//类0（前景）灰度均值
		float class_0_average = 0;
		//类1（背景）灰度均值
		float class_1_average = 0;
		for(int i = 0; i <= k; i++){
			class_0_average += i * (pixelProb[i] / (class_0_prob+0.000001));
		}
		for(int i = k+1; i <= 255; i++){
			class_1_average += i * (pixelProb[i] / (class_1_prob+0.000001));
		}

		//类0（前景）方差
		float class_0_variance = 0;
		//类1（背景）方差
		float class_1_variance = 0;
		for(int i = 0; i <= k; i++){
			class_0_variance += std::pow((i-class_0_average),2)*(pixelProb[i]/(class_0_prob+0.000001));
		}

		for(int i = k+1; i <= 255; i++){
			class_1_variance += std::pow((i-class_1_average),2)*(pixelProb[i]/(class_1_prob+0.000001));
		}

		//类内方差
		float class_within_variance = class_0_prob*std::pow(class_0_variance,2) + class_1_prob*std::pow(class_1_variance,2);
		class_within_variance_vec.push_back(class_within_variance);
		//类间方差
		float class_between_variance = class_0_prob * class_1_prob * std::pow((class_1_average-class_0_average),2);
		class_between_variance_vec.push_back(class_between_variance);
	}

	auto max_value_iter = std::max_element(class_between_variance_vec.begin(), class_between_variance_vec.end());
	int k_value = std::distance(class_between_variance_vec.begin(), max_value_iter);
	std::cout<<"k_value:"<< k_value<<"  max_var:"<<*max_value_iter<<std::endl;

/*	auto min_value_iter = std::min_element(class_within_variance_vec.begin(), class_within_variance_vec.end());
	int kk_value = std::distance(class_within_variance_vec.begin(), min_value_iter);
	std::cout<<"kk_value:"<< kk_value<<"  max_var:"<<*min_value_iter<<std::endl;*/

	for(int i = 0; i < tempImage.rows; i++){
		for(int j = 0; j < tempImage.cols; j++){
			if(tempImage.at<uchar>(i, j) <= k_value)
				tempImage.at<uchar>(i, j) = 0;
			else
				tempImage.at<uchar>(i, j) = 255;
		}
	}

	cv::imshow("lenna", originImage);

	cv::threshold(tempImage_2, tempImage_2, 0, 255, CV_THRESH_OTSU);
	cv::imshow("OpenCV_OSTU", tempImage_2);

	cv::imshow("my_achieve_ostu", tempImage);
	cv::waitKey(0);
}








